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ABSTRACT 

£\\ , We describe a general radiative equilibrium and temperature correction procedure 

for use in Monte Carlo radiation transfer codes with sources of temperature-independent 
opacity, such as astrophysical dust. The technique utilizes the fact that Monte Carlo 
simulations track individual photon packets, so we may easily determine where their 
i-^ ■ energy is absorbed. When a packet is absorbed, it heats a particular cell within the 

envelope, raising its temperature. To enforce radiative equilibrium, the absorbed 



o 



packet is immediately re-emitted. To correct the cell temperature, the frequency of 



the re-emitted packet is chosen so that it corrects the temperature of the spectrum 
previously emitted by the cell. The re-emitted packet then continues being scattered, 
absorbed, and re-emitted until it finally escapes from the envelope. As the simulation 
runs, the envelope heats up, and the emergent spectral energy distribution (SED) 
relaxes to its equilibrium value, without iteration. This implies that the equilibrium 
temperature calculation requires no more computation time than the SED calculation 
of an equivalent pure scattering model with fixed temperature. In addition to avoiding 
iteration, our method conserves energy exactly, because all injected photon packets 
eventually escape. Furthermore, individual packets transport energy across the entire 
system because they are never destroyed. This long-range communication, coupled 
with the lack of iteration, implies that our method does not suffer the convergence 
problems commonly associated with A-iteration. To verify our temperature correction 
procedure, we compare our results to standard benchmark tests, and finally we present 
the results of simulations for two-dimensional axisymmetric density structures. 



Subject headings: radiative transfer — scattering — stars: circumstellar matter - 
ISM: dust, extinction 
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1. Introduction 

There is an ever increasing wealth of observational evidence indicating the non-sphericity 
of almost every type of astronomical object (e.g., extended circumstellar environments, novae 
shells, planetary nebulae, galaxies, and AGNs). To accurately interpret this data, detailed two- 
and three-dimensional radiation transfer techniques are required. With the availability of fast 
workstations, many researchers are turning to Monte Carlo techniques to produce model images 
and spectra for the asymmetric objects they are investigating. In Monte Carlo radiation transfer 
simulations, packets of energy or "photons" are followed as they are scattered and absorbed within 
a prescribed medium. One of the features of this technique is that the locations of the packets are 
known when they are absorbed, so we can determine where their energy is deposited. This energy 
heats the medium, and to conserve radiative equilibrium, the absorbed energy must be reradiated 
at other wavelengths, depending on the opacity sources present. Tracking these photon packets, 
while enforcing radiative equilibrium, permits the calculation of both the temperature structure 
and emergent spectral energy distribution (SED) of the envelope. The ability of Monte Carlo 
techniques to easily follow the transfer of radiation through complex geometries makes them very 
attractive methods for determining the temperature structure within non-spherical environments 
- a task which is very difficult with traditional ray tracing techniques. 

Previous work on this problem for spherical geometries includes the approximate solutions by 
Scoville & Kwan (1976), who ignored scattering, Leung (1976), and diffusion approximations by 
Yorke (1980). The spherically symmetric problem has been solved exactly by Rowan- Robinson 
(1980), Wolfire & Cassinelli (1986), and Ivezic & Elitzur (1997), who used a scaling technique. 
Extensions of the exact solution to two dimensions have been performed by Efstathiou & 
Rowan-Robinson (1990, 1991), while approximate two-dimensional models have been presented by 
Sonnhalter, Preibisch, & Yorke (1995) and Men'shchikov & Henning (1997). 

Radiative equilibrium calculations using the Monte Carlo technique have been presented by 
Lefevre, Bergeat, & Daniel (1982); Lefevre, Daniel, & Bergeat (1983); Wolf, Henning, & Secklum 
(1999); and Lucy (1999). Most of these authors (Lucy being exceptional) use a technique in which 
stellar and envelope photon packets are emitted separately. The number of envelope packets to be 
emitted is determined by the envelope temperature, while the envelope temperature is determined 
by the number of absorbed packets. Consequently these techniques require iteration, usually 
using the absorbed stellar photons to provide an initial guess for the envelope temperature. The 
iteration proceeds until the envelope temperatures converge. Note that the stellar luminosity is 
not automatically conserved during the simulation; only after the temperatures converge is the 
luminosity approximately conserved. 

In contrast, Lucy adopts a strategy in which absorbed photon packets are immediately 
re-emitted, using a frequency distribution set by the current envelope temperature. Although 
the frequency distribution of the reprocessed photons is incorrect (until the temperatures have 
converged), his method automatically enforces local radiative equilibrium and implicitly conserves 
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the stellar luminosity. The insight of Lucy's method is that conservation of the stellar luminosity 
is more important than the spectral energy distribution when calculating the radiative equilibrium 
temperatures. Nonetheless, this method requires iteration. 

The primary problem faced by Lucy's method is the incorrect frequency distribution of the 
re-emitted photons. In this paper we develop an adaptive Monte Carlo technique that corrects 
the frequency distribution of the re-emitted photons. Essentially, our method relaxes to the 
correct frequency and temperature distribution. Furthermore it requires no iteration as long as 
the opacity is independent of temperature. Such is the case for astrophysical dust. In Section 2, 
we describe the temperature correction algorithm. We compare the results of our code with a 
spherically symmetric code in Section 3, and in Section 4 we present results for two dimensional 
axisymmetric density structures. 

2. Monte Carlo Radiative Equilibrium 

We wish to develop a method to calculate the temperature distribution throughout an 
extended dusty environment for use with Monte Carlo simulations of the radiation transfer. The 
radiation transfer technique we employ has been described in detail in other papers: Code & 
Whitney (1995); Whitney & Hartmann (1992, 1993); Wood et al. (1996), so we only summarize 
it here. The basic idea is to divide the luminosity of the radiation source into equal-energy, 
monochromatic "photon packets" that are emitted stochastically by the source. These packets are 
followed to random interaction locations, determined by the optical depth, where they are either 
scattered or absorbed with a probability given by the albedo. If the packet is scattered, a random 
scattering angle is obtained from the scattering phase function (differential cross section). If 
instead the packet is absorbed, its energy is added to the envelope, raising the local temperature. 
To conserve energy and enforce radiative equilibrium, the packet is re-emitted immediately at a 
new frequency determined by the envelope temperature. These re-emitted photons comprise the 
diffuse radiation field. After either scattering or absorption plus reemission, the photon packet 
continues to a new interaction location. This process is repeated until all the packets escape the 
dusty environment, whereupon they are placed into frequency and direction-of-observation bins 
that provide the emergent spectral energy distribution. Since all the injected packets eventually 
escape (either by scattering or absorption followed by reemission), this method implicitly conserves 
total energy. Furthermore it automatically includes the diffuse radiation field when calculating 
both the temperature structure and the emergent spectral energy distribution. 

We now describe in detail how we calculate the temperature structure and SEDs of dusty 
environments illuminated by a radiation source. This radiation can come from any astrophysical 
source, either internal or external, point-like or extended. 
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2.1. Radiative Equilibrium Temperature 

Initially we divide the source luminosity, L, into iV 7 photon packets emitted over a time 
interval At. Each photon packet has the ScimG energy, Ery , so 

£ 7 = LAt/Nry . (1) 

Note that the number of physical photons in each packet is frequency-dependent. 

When the monochromatic photon packet is injected into the envelope, it will be assigned 
a random frequency chosen from the spectral energy distribution of the source. This frequency 
determines the dust absorptive opacity, k u , and scattering opacity, a u (both per mass), as well as 
the scattering parameters for the ensuing random walk of the packet through the envelope. The 
envelope is divided into spatial grid cells with volume Vi, where i is the cell index. As we inject 
source photon packets, we maintain a running total, iV», of how many packets are absorbed in 
each grid cell. Whenever a packet is absorbed in a grid cell, we deposit its energy in the cell and 
recalculate the cell's temperature. The total energy absorbed in the cell 



We assume that the dust particles are in local thermodynamic equilibrium (LTE), and for 
simplicity we adopt a single temperature for the dust grains. Note that although we use dust for 
the continuous opacity source, we could replace the dust by any continuous LTE opacity source 
that is independent of temperature. In radiative equilibrium, the absorbed energy, Ef hs , must 
be reradiated. The thermal emissivity of the dust j u = K u pB u (T), where B U (T) is the Planck 
function at temperature T, so the emitted energy is 

Ef m = 4vrAt J dVi J P K u B u (T)du 

= 4vrAt J n P (T)B(T)pdVi , (3) 

where up = J k v B v dvj B is the Planck mean opacity, and B = aT 4 /it is the frequency integrated 
Planck function. If we adopt a temperature that is constant throughout the grid cell, Tj, then 

Ef™ = AirAtKpiT^BiT^mi , (4) 

where mj is the mass of the cell. 

Equating the absorbed (||) and emitted (|j) energies, we find that after absorbing iVj packets, 
the dust temperature is given by 

T 4 = NjL 

a 1 4iV 7Kp (T i )m i ' 1 1 

Note that the Planck mean opacity, Kp, is a function of temperature, so equation (0) is actually 
an implicit equation for the temperature, which must be solved every time a packet is absorbed. 
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Since this equation is solved so many times, an efficient algorithm is desirable. Fortunately Kp 
is a slowly varying function of temperature, which implies a simple iterative algorithm may be 
used to solve equation (||). To do so efficiently, we pre-tabulate the Planck mean opacities for a 
large range of temperatures and evaluate Kp(Tj) by interpolation, using the temperature from the 
previous guess. After a few steps, we have the solution for T. Note that because the dust opacity 
is temperature-independent, the product KpaTf, which is proportional to J K u B u dv, increases 
monotonically with temperature. Consequently T always increases when the cell absorbs an 
additional packet. 



2.2. Temperature Correction 

Now that we know the temperature after absorbing an additional packet within the cell, we 
must reradiate this energy so that the heating always balances the cooling. Prior to absorbing this 
packet, the cell previously has emitted packets that carried away an energy corresponding to the 
cell's previous emissivity j' v = K v B v (Ti — AT), where AT is the temperature increase arising from 
the packet absorption. Note that these previous packets were emitted with an incorrect frequency 
distribution corresponding to the previous temperature, Tj — AT. The total energy that should 
be radiated now corresponds to j v at the new temperature, T. Thus the additional energy to be 
carried away is given by 

Aj„ = ju ~ Ju = [Bu(Ti) ~ B u (Ti - AT)] , (6) 

which is the shaded area shown in Figure |l|. As long as the packet energy T 7 is not too large 
(this may be assured by choosing a large enough number of photon packets, 2V~, to use for the 
simulation), the temperature change AT is small, so the temperature correction spectrum 

Aj v « ^AT^ . (7) 

Note that Aj u is everywhere positive because AT > 0, and the Planck function is a monotonically 
increasing function of temperature. Therefore to correct the previously emitted spectrum, we 
immediately re-emit the packet (to conserve energy), and we choose its frequency using the shape 
of Aj u . This procedure statistically reproduces Aj u for the distribution of the re-emitted packets. 
Normalizing this distribution, we find the temperature correction probability distribution 



dPi k v ( dB L 



dv K \ dT ) t=t 



(8) 



where dPi/dv is the probability of re-emitting the packet between frequencies v and v + dv, and 
the normalization constant K = K u (dB u /dT) dv. 

Now that the packet's frequency has changed, we change the opacity and scattering parameters 
accordingly and continue with scattering, absorption, temperature correction, and re-emission 



1 



10 100 



X ((a.m) 

Fig. 1. — Temperature Correction Frequency Distribution. Shown are the dust emissivities, 
ju = HuB u (T), prior to and after the absorption of a single photon packet. The spectrum of the 
previously emitted packets is given by the emissivity at the old cell temperature (bottom curve). To 
correct the spectrum from the old temperature to the new temperature (upper curve), the photon 
packet should be re-emitted using the difference spectrum (shaded area). 
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until all the photon packets finally escape from the system. In principle we could also account for 
back-warming of the source. Whenever a packet hits the source, the source must reradiate this 
new energy. This will change the temperature of the source, and the new source photons can be 
emitted using a difference spectrum similar to equation (||). 

When we begin our calculation, no packets have been absorbed, so the initial temperature is 
zero throughout the envelope. This means that the initial temperature change is not small as is 
required by equation ([?]). One could use equation (g) to re-emit the first packet that is absorbed; 
however, this is not necessary. The number of packets producing this initial temperature change 
is small; it is of order the number of spatial grid cells. Furthermore, these packets generally are 
re-emitted at such long wavelengths that they are not observable. Consequently the error arising 
from using equation (|7]) to re-emit every packet is too small to be of importance. 

As the simulation runs, the envelope starts at zero temperature. It then heats up, and the 
radiation field "relaxes" to its final spectral shape. The temperature correction procedure is simply 
a way of re-ordering the calculation (which frequencies are being used at a given moment) so that 
in the end, all the frequencies have been properly sampled. Consequently, after all the stellar 
photon packets have been transported, the final envelope temperature is the correct radiative 
equilibrium temperature, and the emergent spectral energy distribution has the correct frequency 
shape. Furthermore, the photon re-emission automatically accounts for the diffuse envelope 
emission. Note that energy is necessarily conserved, and there is no time-consuming iteration in our 
scheme; calculating the radiative equilibrium temperature requires no more computational time 
than an equivalent pure scattering calculation in which the temperature structure is held fixed. 
Similarly, there is no issue of convergence in our method. Unlike A-iteration, the photon packets 
carry energy over large distances throughout the envelope because they are never destroyed, and 
of course there is no iteration at all. After running N-y packets, we have the final answer. The 
simulation does not continue running until some convergence criterion is met, and the only source 
of error in the calculation is the statistical sampling error inherent in Monte Carlo simulations. 



3. Benchmark Verification 

To validate our method for determining the radiative equilibrium temperatures and emergent 
fluxes, we compared our results against a set of benchmark calculations recently developed by 
Ivezic et al. (1997) for testing spherically symmetric dust radiative equilibrium codes. The 
parameters listed by Ivezic et al. enable us to exactly reproduce the same set of physical conditions 
(i.e, input spectrum, dust destruction radius, optical depth, opacity frequency distribution, and 
radial density structure). 

For all benchmark tests, Ivezic et al. used a point source star, radiating with a black body 
spectrum whose temperature T* = 2500K. The dust density distribution was a power law with 
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radius, 
where 

= Tx fl/^max = 0), 

P ° (Kx + a x )(l-R dnst /R m&x )\l/R dnst = 2). 1 > 

The inner radius of the envelope is the dust destruction radius, i?dust> the outer radius is -Rmax, 
and the total radial optical depth is t\, specified at A = lum. The dust absorptive opacity, k u , 
and scattering opacity, a u , were taken to be 

K v _ 05 /l (A<l/im), 



+ cr)i Mm I (lum/A) (A>lum), 

a v { 1 (A < 1 urn), , s 

°- 5 i ^.._,xx4 /x . i .~~\ (11) 



(k + c)i^ m I (lum/A) 4 (A>lum). 

Since the total optical depth at lum is independently specified, these opacities have been 
normalized to that at lum for convenience. The wavelength-dependent scattering albedo is given 
by a = ov/O^ + ov), and the scattering was assumed to be isotropic (note that in dust simulations, 
we would normally use a non-isotropic phase function for the scattering). 

In principle, the inner radius of the dust shell, i?dust> is determined by the dust condensation 
temperature, chosen by Ivezic et al. to be T cond = 800 K. However, we are only testing the 
temperature correction procedure, so we have not implemented a scheme to solve self consistently 
for the dust destruction radius. The values for .Rdust were calculated instead using Ivezic et 
al.'s eq. (4) and data from their Table 1. Finally, the outer radius of the dust shell was set to 
be -R max = 10 3 i?dust- The parameters describing the various test simulations are summarized in 
Table @. 

To begin the simulation, we release stellar photon packets with a black body frequency 
distribution, given by the normalized Planck function 

, , s 15 x 3 , . 

7r 4 e x — 1 

where x = hh>/kT+. A particularly simple method for sampling the black body distribution is given 
by Carter & Cashwell (1975). Since this reference is somewhat obscure and difficult to obtain, we 
summarize the method here. First, generate a uniform random number, £o> in the range to 1, 
and determine the minimum value of I, l m \ n , that satisfies the condition 

i=i 

Next obtain four additional uniform random numbers (in the range to 1), £i, £2, £3; and £4. 
Finally, the packet frequency is given by 



x = - In (£l6£3£4)/Jmin 



(14) 
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Table 1. Spherical Models 



p 




-Rdust/-R* 





1 


8.44 





10 


8.46 





100 


8.60 


2 


1 


9.11 


2 


10 


11.37 


2 


100 


17.67 



Table 2. Ellipsoidal Models 



Peq/ Ppole 


-Rdust/ R* 


eq 

V 


pole 
T V 


1000 


10 


200 


0.2 


1000 


10 


20 


0.02 


10 


10 


200 


20 


10 


10 


20 


2 
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After emitting these packets from the star, we track them through the envelope. To determine 
the envelope temperature, we must count how many packets are absorbed in each grid cell. Since 
the envelope is spherically symmetric, we employ a set of spherical shells for our grid. To obtain 
the best Poisson error statistics, we should ideally construct the grid positions so that equal 
numbers of packets are absorbed in each cell. Since the probability of absorbing a photon packet 
is proportional to the optical depth, we choose equal radial optical depth grid locations, 

iAr + 1 (p = 0), 

(15) 



i^dust [Nf/(N f -i) (p = 2), 

where N r = 200 is the total number of cells we used, Ar = (-R m ax/-Rdust — l)/N r , and 

N f = N r /(1 - R dust /R max ). Integrating the density, eq. (|9|), over the cell volume to obtain the 

mass, we find from eq. (0) that the temperature in each grid cell is given by 

! A f »A f r(i?*/-Rd U st) 2 [(^lMm+^l M m)/Kp(^)] , _ ^ 

4A r 7 T lMm Ki 2 -i+l/3)Ar^+(2i-l)Ar+lj W U J' , s 

^(Af f -i)(Af / -i+l)( J R,/fidu St ) 2 (l-«d ust /fi m ax) / _ 9 x {L0) 

4Af 7 r lAlm 7V / ^ F (T i )/(Ki flm +CTi Mm )l ^ P » 

where we have used L = AirR^aT^ for the stellar luminosity. 

We then proceed with the radiation transfer, temperature calculation, and reemission as 
described in Section 2 until all packets exit the envelope. When the packets escape, they are 
placed into N u = 1024 uniform frequency bins, 

v k = ^maxTT , (17) 

where /iz^max = 16kT*. The width of each bin Av = Umax/Ny. Since the envelope is spherically 
symmetric, the observed flux F Uk = N^E^ / Aird 2 At/S.v , where is the number of packets in the 
k th frequency bin, and d is the observer's distance from the star. Normalizing to the total flux, 
F = L/Aird 2 , the SED is given by 

t)r«- 1/2) w,- < 18) 

The factor (k — 1/2) arises from using the frequency at the center of the bin. 



In Figures [2a| and 2b we show the results of our simulation compared with the output of one 
of the codes tested by Ivezic et al. This code, called DUSTY, is publicly available and is described 
in Ivezic & Elitzur (1997). We see that our Monte Carlo calculations reproduce both the correct 
temperature structure and SED. Note the large errors at the longest and shortest wavelengths 
in the Monte Carlo calculations. At these wavelengths, the number of emerging packets is 
small, resulting in large errors (the relative error in each flux bin equals lfy/N^, owing to the 
Poisson sampling statistics inherent in Monte Carlo simulations). The excellent agreement of the 
comparisons shown in Figures 2a and 2b (the differences are smaller than the numerical error 
of the DUSTY calculations) demonstrates the validity of our temperature correction procedure 
described in Section 2. 
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X (|im) 

Fig. 2a. — Spherical Model SEDs. Shown are the comparisons of our spherically symmetric model 
fluxes (histogram) with those from DUSTY (solid curves) . The top panel compares the 1/r 2 density 
models, while the bottom panel compares the constant density models. The optical depths at 1/im 
are as indicated. The incident stellar spectrum is shown by the dashed curve. 
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Fig. 2b. — Spherical Model Temperatures. Shown are the comparisons of our spherically symmetric 
model temperatures (histogram) with those from DUSTY (solid curves). The top panel compares 
the 1/r 2 density models, while the bottom panel compares the constant density models. The optical 
depths at 1/xm are as indicated. 
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Now that we have verified our our basic radiative equilibrium algorithm, we can proceed to 
investigate the temperature structure and SEDs of other geometries. Owing to the inherently 
three-dimension nature of Monte Carlo simulations (even our 1-D spherically symmetric code 
internally tracks the photon packets in three dimensions), our code is readily modified for 
arbitrary geometries. We now show the results of an application to axisymmetric circumstellar 
environments. 



4. Axisymmetric 2-D Calculations 

For the purpose of this illustrative simulation, we adopt a stellar blackbody with = 3500 K, 
envelope inner radius iZdust = 10i?*, and a simple ellipsoidal parameterization for the circumstellar 
density. The isodensity contours are elliptical with a/b being the ratio of the semi-major to 
semi-minor axis. The density is given by 

(RdustM 2 , . 

1 + j z cos z 

where po is given by eq. ([n]; p = 2 case), 6 is the polar angle, and the "flattening factor" 



/ = y/(a/b)* - 1 • (20) 
Note that the equatorial to polar density ratio at a given radius is p eq / p po \ e = (a/b) 2 . 

As before, we divide the circumstellar environment into cells with N r = 200 radial and 
= 20 latitudinal grid points. Note that the envelope is symmetric about the equator, so 

we combine the cells below the equator with their counterparts above the equator. This is 

automatically accomplished by using p,j = cos(9j) for the latitudinal grid point coordinate. 

Spacing the grid so that the radial and latitudinal optical depths are the same for each cell, we 

find 

N f 

N = i 1 (21) 

.yi + (a/&) 2 tan 2 (7rj/2Ay 

With these cell coordinates, equation (ph for the temperature in cell becomes 

T 4 _ T 4 N id (N f -i)(N f -i + 1)(R*/R dust ) 2 (l - R dnst /R max )f 
" MA^ 7 r^ / [ Kp (r^O/(^ + ^)][tan-HM-i)-tan-i(M)] ' 1 ' 

where Nij is the number of packets absorbed in the cell, and we have chosen A = 5500A (l/-band) 
for our equatorial radial optical depth parameter, ri 



cq 

v ■ 



For the dust opacity, we adopt a standard MRN interstellar grain mixture (Mathis, Rumpl, & 
Nordsieck 1977), using optical constants from Draine & Lee (1984). Figure [|| shows the opacity and 
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X (urn) 

Fig. 3. — Dust Opacity. The normalized absorptive (dashed line), scattering (dotted line), and 
total (solid line) opacities for the 2-D simulations are shown in the top panel. The bottom panel 
shows the the corresponding albedo. 
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albedo in graphical form. The data for this figure is available in tabular formQ from the DUSTY 
Web site, http: / / www.pa.uky.edu/~moshe/dusty . Note the prominent silicate absorption features 



at 10pm and 18/iin. For the current demonstration, we assume that the scattering is isotropic, 
but we can easily accommodate any phase function, analytic or tabulated when calculating the 
emergent SED. 

Unlike the spherically symmetric benchmarks, the SED now depends on viewing angle, so we 
must bin the escaping packets both in direction and frequency. Since the envelope is axisymmetric, 
the observed flux only depends on the inclination angle, i, of the envelope symmetry axis (i.e., 
the stellar rotation axis). To obtain approximately equal numbers of escaping packets, we 
choose A^i nc = 10 direction bins, with equal solid angles, Af2 = 47r/iVi nc . The escaping directions 
(inclination angles), ij, for these bins are given by 

Mf C = l/Ninc , (23) 

where pf* = cos(^). In addition to the direction bins, we use the frequency bins given by eq. (|l7|). 

Now that we have specified both the direction and frequency bins, the observed flux is 
Fy kl = N^iEy/ AQcPAtAv, where N^i is the number of escaping packets in the (k,l) bin. 
Normalizing to the bolometric flux F gives the emergent SED, 

( vF v \ _ {k — l/2)N k jNinc (24) 



V F J kjl N 7 

We now choose to investigate the SEDs produced by two density structures with different 
degrees of flattening. The first has a density ratio p eq / p po i e = 1000 to represent a disk-like 
structure. The second has a density ratio p eq / p vo \ e = 10, which is mildly oblate, representative 
of an infalling protostellar envelope. For each density structure, we perform optically thick (in 
the mid-IR) and optically thin calculations. The optically thick calculations have an equatorial 
1^-band optical depth Ty q = 200, while the optically thin calculations have Ty q = 20. Table [2] 
summarizes these density structures. For comparison, we also have performed calculations for 
p = 2 spherically symmetric models containing the same total mass as our disk and envelope 
densities. Keeping the same inner and outer radii, the optical depth for the equivalent spherical 
model is 

4 P = rf-^-l . (25) 



4.1. Disk Model 



Figure 4a shows the incident stellar spectrum and the emergent SED as a function of 
viewing angle for the disk-like model, as well as the result of the equivalent spherically symmetric 



1 ftp:/ /gradj .pa.uky.edu / dusty/distribution / ism-stnd.dat 
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1 10 100 

X (jim) 



Fig. 4a. — Disk Model SEDs. The normalized emergent fluxes (solid lines) are shown as a function 
of viewing angle (10 inclinations, evenly spaced in cosi) for the disk-like circumstellar density. 
The lowest curve corresponds to an almost edge-on view (cosi = 0.05), while the highest curve 
corresponds to an almost pole-on view (cosi = 0.95). The dotted curves are for a spherically 
symmetric simulation, having the same total circumstellar mass. The incident stellar spectrum is 
shown by the dashed curve. 



-17- 




1 10 100 1000 



r/R H , 

dust 

Fig. 4b. — Disk Model Temperatures. The temperatures (solid lines) are shown as a function of 
polar angle (20 angles, spaced as indicated in the text). The lowest curve is closest to the equatorial 
plane and the highest curve is nearest to the polar axis. The dashed curve is the temperature for 
a spherically symmetric simulation, having the same total circumstellar mass, and the light dotted 
curve is a simple power law T oc r~ 0A for reference. 
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calculation. For both disk simulations, when viewing the system pole-on, we are looking through 
optically thin circumstellar dust (see Table |2|) and can see the star at optical wavelengths. In the 
IR, there is an excess arising from the circumstellar disk, which reprocesses the stellar radiation. 
Note that, for pole-on viewing, we see the silicate features in emission, since the disk is optically 
thin in the vertical direction. As we go to higher viewing angles, the optical depth to the central 
star increases, and as a result, the star becomes more extincted in the optical region. Note that 
at almost edge-on viewing, a "shoulder" appears around 1/im for the optically thick case. This 
arises due to the dominant effects of scattering the stellar radiation at these wavelengths. These 
scattering shoulders are also present in the axisymmetric calculations presented by Efstathiou & 
Rowan-Robinson (1991), Sonnhalter et al. (1995), Men'schikov & Henning (1997), and D'Alessio 
et al. (1999). At wavelengths longer than about 1/xm, the albedo begins to drop rapidly (see 
Fig. ||), and the disk thermal emission begins to dominate, so the shoulder terminates. Beyond 
30/xm, the envelope becomes optically thin, so the spectrum is independent of inclination and is 
dominated by the dust emission. 

The two dimensional temperature structure for the disk-like models is displayed in Figure 4b. 
We see that at the inner edge of the envelope (i?dust) there is little variation of the temperature 
with latitude, while at large radii there is a clear latitudinal temperature gradient, with the dust 
in the denser equatorial regions being cooler than dust at high latitudes. 

In the polar region, the material is optically thin to the stellar radiation, so it heats up to 
the optically thin radiative equilibrium temperature. This temperature has a power law behavior, 



T oc r 0,4 for dust opacity (ft tx A l ). As can be seen in Figure 4b, the polar temperature does 



indeed have a power law decrease with a slope of approximately —0.4. 

In contrast, the disk only displays this power law behavior at large radii. At the inner edge of 
the disk, the disk sees the same mean (stellar) intensity as is present in the polar region (J = WB, 
where W = 0.5{1 — [1 — (i?/r) 2 ] 1//2 } is the dilution factor). Consequently, the inner edge of the 
disk heats up to the optically thin radiative equilibrium temperature — the same as the polar 
temperature. However at larger radii, the opaque material in the inner regions of the disk shields 
the outer regions from direct heating by the stellar radiation. This shielding reduces the mean 
intensity (J < WB). As a result, the outer disk is only heated to a fraction of the optically thin 
radiative equilibrium temperature. Thus the equatorial region is cooler than the polar region. 
Eventually at large enough radii, the disk becomes optically thin to the heating radiation. At 
that point, it sees a radially streaming radiation field from an effective photosphere that has a 
much lower temperature than the star. From that point outward, the disk temperature displays 
an optically thin power law decrease with a slope that parallels the polar temperature. 

The spherically symmetric calculation overestimates both the emergent flux and the disk 
temperature. In the optically thin limit, the IR continuum is proportional to the mass of the 
circumstellar material, so one would expect that a spherically equivalent mass would reproduce the 
long wavelength spectrum. Recall however that we can see the star at pole-on viewing angles. This 
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implies that the disk does not reprocess the entire bolometric luminosity of the star. Consequently, 
the IR excess is less than that in the spherical model. Note that one can nonetheless reproduce the 
edge-on SED using a spherical model if we allow the density power law to depart from 1/r 2 and 
change the size of the spherical envelope. This also was noted by Sonnhalter et al. (1997), who 
found that they could fit the IR continuum of their disk models by changing the radial dependence 
of the circumstellar density and the outer radius for the same total envelope mass. 

4.2. Envelope Model 

The spectral energy distribution and temperature structure for the oblate envelope are shown 
in Figures [5a] and [5t]. The optically thin envelope displays significant extinction of the star at all 
viewing angles. Close to edge-on, scattering shoulders appear around l//m that are similar to those 
seen in the disk model (Fig. ^a|). Finally, because the model is optically thin in the mid-IR, the 
silicate feature is always in emission, and the SED is independent of viewing angle for wavelengths 
longer than a few microns. 

For the denser envelope, the star is extremely faint at optical wavelengths. Along with the 
increased extinction in the optical, the scattering shoulders are less prominent, due to the thermal 
emission by the envelope. In the mid-IR, the envelope is optically thick edge-on and optically thin 
pole-on. Consequently, the silicate features go from absorption to emission as the viewing angle 
changes from edge-on to pole-on. The general shape of the SEDs, which now peak in the mid-IR 
for all inclinations, are reminiscent of spectra from embedded (Class I) T Tauri stars , which are 
commonly modeled using flattened axisymmetric dusty envelopes (e.g., Adams, Lada, & Shu 
1987; Kenyon, Calvet, & Hartmann 1993; Men'shchikov & Henning 1997; D'Alessio, Calvet, & 
Hartmann 1997). Of these T Tauri simulations, only Efstathiou & Rowan-Robinson (1991) have 
performed an exact calculation for the SED and circumstellar temperature for the Terebey, Shu, 
& Cassen (1984) collapse model. 

The temperature structure for our envelope models (Fig. ^) is qualitatively similar to the 
temperature structure of the disk models (Fig. ^t]). The temperature at the inner edge of the 
envelope is independent of latitude, while at larger radii the equatorial regions are cooler than 
the polar regions. The primary difference is that the latitudinal temperature gradient is not as 
extreme. 

Finally, we note that the equivalent spherically symmetric model better reproduces the far 
IR spectrum of the envelope models than the disk models. This is because the envelope models 
reprocess the entire bolometric luminosity, while in disk models, some of the stellar luminosity 
escapes through the polar region. 
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Fig. 5b. — Envelope Model Temperatures. Same as Fig. 4b, but for the envelope density- 
distribution. 
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5. Discussion 

We have developed a temperature correction procedure for use in Monte Carlo radiation 
transfer codes. We have tested our method against other spherically symmetric benchmark codes 
and successfully matched their results. After verifying our method, we applied it to obtain sample 
temperature distributions and SEDs for 2-D axisymmetric disk-like models and mildly oblate 
envelopes. These simulations illustrate the important role envelope geometry can play when 
interpreting the SEDs of embedded sources. 

The primary limitation of our temperature correction procedure is that it only applies if 
the opacity is independent of temperature. This is not true for free-free opacity or hydrogen 
bound-free opacity, but it is true for dust opacity, which is the dominant opacity source in many 
astrophysical situations. The reason our method is limited to temperature-independent opacities 
is that two associated problems occur if the opacity varies when the cell's temperature changes. 
The first is that the cell will have absorbed either too many or too few of the previous packets 
passing through the cell. The second is the associated change in the interaction locations of the 
previous packets, which implies that the paths of the previous photon packets should have been 
different. These problems do not occur if the opacity is independent of temperature. 

Lucy (1999) has proposed a slightly different method to calculate the equilibrium temperature. 
Instead of sampling photon absorption, he directly samples the photon density (equivalent to the 
mean intensity) by summing the path length of all packets that pass through a cell. Typically 
more packets pass through a cell than are absorbed in the cell, so this method potentially produces 
a more accurate measurement of the temperature for a given total number of packets, especially 
when the envelope is very optically thin. The disadvantage of his method is that it requires 
iteration to determine the envelope temperature. In principle one could partially combine both 
methods. First, run our simulation, adding pathlength counters to each cell. After running all 
packets, use the pathlength information to calculate a final temperature. This will provide a more 
accurate temperature that can be used to calculate the source function. After obtaining the source 
function for each cell, it is a simple matter to integrate the transfer equation to obtain the SED. 

Another limitation of the Monte Carlo method is that it is not well suited to envelopes with 
very high optical depths (r > 100 - 1000), unless there is an escape channel for the photons. 
For example, geometrically thin disks can be extremely optically thick in the radial direction 
but optically thin in the polar direction, providing an escape channel for the photons. Similarly, 
dense envelopes can be very optically thick to the illumination source, but optically thin to the 
reprocessed radiation, which is a another escape channel. In the event no such escape channels 
exist, one must turn to other methods. We are currently investigating how to couple Monte Carlo 
simulation in the optically thinner regions with other methods in the optically thick interior. 

The focus of the present paper has been on the development and implementation of the 
temperature correction procedure. For this reason, we have made several simplifying assumptions 
regarding the circumstellar opacity and geometry. For example, we have assumed a single 
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temperature for all dust grains regardless of their size and composition. This differs from other 
investigations where different types of grains can have different temperatures at the same spatial 
location (e.g., the spherically symmetric radiation transfer code developed by Wolfire & Cassinelli 
1986). Similarly, we have not implemented a procedure to solve for the location of the dust 
destruction radius, which will be different for grains of differing size and composition. These issues 
are currently under investigation. 

With the speed of today's computers, Monte Carlo radiation transfer simulations can be 
performed in a reasonably short time. The 2-D simulations presented in this paper employed 10 8 
packets, requiring about two hours of CPU time; the spherically symmetric cases used 10 6 packets, 
requiring about one minute. So for continuum transfer, Monte Carlo simulation is proving to be a 
very powerful technique for investigating arbitrary density structures and illuminations. 

We would like to thank Barbara Whitney and Mario Magalhaes for many discussions relating 
to this work. An anonymous referee also provided thoughtful comments which helped improved the 
clarity of our presentation. This work has been funded by NASA grants NAG5-3248, NAG5-6039, 
and NSF grant AST-9819928. 
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